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ABSTRACT 

We propose and implement a fast, universally applicable method for extracting the angular power 
spectrum Cn from CMB temperature maps by first estimating the correlation function ^ {9) . Our proce- 
dure recovers the C^'s using N'^ (but potentially A^logiV), operations, where N is the number of pixels. 
This is in contrast with standard maximum likelihood techniques which require operations. Our 
method makes no special assumptions about the map, unlike present fast techniques which rely on sym- 
metries of the underlying noise matrix, sky coverage, scanning strategy, and geometry. This enables for 
the first time the analysis of megapixel maps without symmetries. The key element of our technique is 
the accurate multipole decomposition of £^{6). The Cg error bars and cross-correlations are found by a 
Monte-Carlo approach. We applied our technique to a large number of simulated maps with Boomerang 
sky coverage in 81000 pixels. We used a diagonal noise matrix, with approximately the same amplitude 
as Boomerang. These studies demonstrate that our technique provides an unbiased estimator of the C^'s. 
Even though our method is approximate, the error bars obtained are nearly optimal, and converged only 
after few tens of Monte-Carlo realizations. Our method is directly applicable for the non-diagonal noise 
matrix. This, and other generalizations, such as minimum variance weighting schemes, polarization, and 
higher order statistics are also discussed. 



1. INTRODUCTION 

Future missions of measuring the Cosmic Microwave 
Background (CMB) fluctuations will revolutionize our 
knowledge of cosmology. They will either confirm or re- 
fute the basic Big Bang paradigm, will reveal the values of 
most cosmological parameters within a few percent (e.g., 
Spergel 1994, Knox 1995, Hinshaw, Bennett, & Kogut 
1995, Jungman et al. 1996, Zaldarriaga, Spergel, & Sel- 
jak 1997, Bond, Efstathiou, & Tegmark 1997). The large 
number of pixels contained in current and future experi- 
ments enable these exciting developments, but at the same 
time, present unprecedented challenges the data analysis. 
The mainstream maximum likelihood techniques are al- 
ready pushed to their limits by the largest existing CMB 
maps, but they are clearly inadequate for future megapixel 
surveys. Therefore the most important near term task for 
CMB research is to find techniques which could perform 
the required analyses with realistic resources, thus fulfill 
the promise of the high precision experiments. 

This Letter proposes a fast, universally applicable 
method for the estimation of C/'s from any CMB pixel 
map, based on estimating correlation functions. The re- 
covered errorbars are at most 10% larger than the theoret- 
ically smallest possible ones imposed by cosmic variance 
and noise. Optimal methods could further decrease these 
errorbars only slightly, and at an unrealistic cost; thus 
they represent a case of diminishing returns. Our pro- 
cedure scales as N"^, and more sophisticated algorithms 
will improve this to TVlogiV. This essentially solves the 
problem for two major upcoming space missions, MAP 
(Microwave Anisotropy Probe) and the Planck, and at 
the same time allows more detailed analyses of current 
and future balloon-borne and ground based experiments. 



The short analysis turn around time enables Monte-Carlo 
(MC) studies of systematics, foregrounds, errors, underly- 
ing models etc. These are essential to the interpretation 
of the data but inaccessible to present methods. 

The standard maximum likelihood methods for analyses 
were developed and tested on COBE measurements (e.g., 
Gorski 1994, Gorski et al. 1994, 1996, Bond 1995, Tegmark 
& Bunn 1995, Hinshaw et al. 1996, Tegmark 1996, Bunn & 
White 1996, Bond, Jaffe, & Knox 1998, hereafter BJK98, 
2000) which have only about N ~ 1000 pixels. Balloon- 
borne experiments, such as Boomerang (de Bernardis et al. 
2000), Maxima (Hanany et al. 2000), and Tophat (Martin 
et al. 1996), and ground based measurements e.g., TOCO 
(Miller et al. 1999) and Viper (Peterson et al. 2000) with 
up to V ~ 10^ are already pushing present day supercom- 
puter technology. The future missions MAP and Planck 
with N = 10^ — lO'^ are estimated to require up to mil- 
lions of years (Borrill 1999, Bond, Crittenden, Jaffe, & 
Knox 2000) for one iteration of the quadratic maximum 
likelihood estimator for C/'s. The disk storage and mem- 
ory usage of these algorithms are prohibitive as well. 

To date fast algorithms were presented under two fairly 
restrictive assumptions; i) the noise is both temporally 
uncorrelated and spatially axially symmetric ii) the fore- 
grounds can be exactly removed both from the map and 
the correlation matrix (Oh, Spergel, & Hinshaw 1999, 
Wandelt, Hivon & Gorski 1998, 2000). Our approach is 
significantly different from these, since it does not assume 
any symmetries. 

2. DESCRIPTION OF THE METHOD 

The basic steps of our recipe to extract C/ 's from a large 
CMB map are conceptually simple: first we measure the 
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Correlation function (MC simulation) 




Fig. 1. — The left panel displays a typical high resolution measurement of the two-point correlation function along with its Gaussian 
resampling at Legendre roots (large dots). The singular feature at around 1 degree (cos 9 ~ 0.9998) condenses most of the information about 
the Ci's. The right panel present a few examples of the inversion process. Individual Cj's are calculated from the correlation function with 
the method described in the paper. 



two-point correlation function in high resolution bins, next 
we smooth it with Gaussian kernels centered on the roots 
of Legendre polynomials, finally, we integrate to obtain 
the C/'s. Each step is fine-tuned in order to arrive at a 
fast method, which is as precise as possible. 

Let us denote the temperature fluctuations at a sky vec- 
tor q, a unit vector pointing to a pixel on the sky, with 
T{q). In isotropic universes the two-point correlation func- 
tion is a function only of the angle between the two vectors 
and can expanded into a Legendre series, 

where cos0 — qi-q2, the dot product of the two unit vec- 
tors, Pi{x) is the ^-th Legendre polynomial, and the Ci 
coefficients realize the angular power spectrum of fluctu- 
ations. If the CMB anisotropy is Gaussian, which is ex- 
pected to be an excellent approximation, the correlation 
function, or the Cgs yield full statistical description. 

In reality each pixel value, A^, contains contributions 
from the CMB and noise; the latter is also assumed to 
be Gaussian with a correlation matrix (the noise matrix) 
Nij, determined during map-making. The full pixel-pixel 
correlation matrix is Cij = + Nij . We adopted the edge 
and noise corrected estimator of Szapudi & Szalay (1998), 

|(cos0) = ^/.,(A,A,-7V,,), (2) 

where /y — unless the pair of pixels belong to a partic- 
ular bin in cos0, and '^^j fij — 1. The above estimator 

is unbiased, i.e. (^(cos0)) = ^(cos0). A complex pair 
weighting might be required in surveys with uneven cov- 
erage, and to optimize the performance of the estimator 
(see discussions in §4). For this Letter we used a straight- 
forward N'^ program to estimate the above quantity in a 
large number of bins linear in 9 massively oversampling the 
pixel separation. The results are robust against the exact 
number of bins used, for the measurements presented here 
we have placed 300, 000 bins in the range of — tt. 



The raw correlation function is then resampled with a 
Gaussian filter with additional weighting proportional to 
the number of pairs in each bin. The centers of the filters 
are determined by the roots of the Legendre polynomial 
^fmax of the highest imax to be measured; their widths is 
about half the pixel size. This procedure allows accurate 
centering of the filters, suppresses any high frequency sub- 
beam measurement noise due to the massive oversampling. 
The Gaussian convolution translates into a simple multi- 
plication in ^-space, thus easy to correct for. We verified 
that, after this correction, the final result is robust against 
varying the size of the filter; we have used 4'. For small 
surveys, a smooth cut-off of the correlation function is nec- 
essary on large angles to suppress noise from edge effects. 
We checked that the Cgs are robust against changing the 
scale and sharpness of the cut-off (up to a corresponding 
smoothing in £-space); we have used exponential tapering 
above 18°. The left panel of Figure |l| shows an example 
of a measured correlation functions along with its resam- 
pled version (large dots). The C^'s are determined from 
the resampled correlation function by Legendre-Gauss in- 
tegration (e.g.. Press et al. 1992). Since this is exact up 
to 2trnax + 1, the Cf's are as accurately recovered as the 
correlation function. The right panel of Figure |l| shows ex- 
amples of 's inverted from two-point correlation function 
measurements in simulations. 

The full correlation matrix of the above estimator can be 
calculated analytically as well. Under Gaussian assump- 
tion, the cross-correlation between two bins denoted with 
1 and 2, respectively is {5^18^2) = J2^jkl fljfkMiki^ji + 
Nji) + k ^ I)]. The straightforward calculation of the 
above expression takes N* operations, thus infeasible for 
a megapixel survey. While this scaling is likely to be im- 
proved to N{logN)^ with the advanced algorithms in the 
near future, in this Letter we propose an entirely different 
approach: high volume MC simulations. Certain system- 
atics can only be taken into account this way, and the 
speed of our method makes such computations possible. 

3. SIMULATIONS 
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Fig. 2. — The left panel displays the main results of this Letter: Ce's calculated in 1298 Boomerang-like simulations and then rebinned 
into flat Cf-bands with width of 50. The small poixits show the individual measurements with errorbars representing the standard deviations 
in each band. Theoretical errorbars of Equation H are displayed shifted to the right for clarity. The comparison shows that our method is 
unbiased, and has nearly optimal variance. The arrows point to the effective beam and pixel scales. The right panel shows the distribution 
of measurements in selected bands. Shifted lognormal (dash) and Gaussian curves are fitted to the histogram. Except for the £ = 50 band, 
the two are indistinguishable. Horizontal lines display optimal FWHM dispersion for comparison. 



To illustrate the method, we generated 1298 MC real- 
izations of CMB maps with a standard CDM power spec- 
trum and typical instrumental noise. The sky coverage 
used to create the maps 1000 deg^), the noise level 
45 /iK per pixel), and the Gaussian beam size (10' 
FWHM) were comparable to the BOOMERanG experi- 
ment (de Bernardis et al. 2000). Full sky maps in ~ T 
pixels were generated by the HEALPix (Gorski et al. 1998) 
synf ast routine. Then maps were reduced to the desired 
coverage, and uniform white noise has been added. 

The angular power spectra have then been extracted, 
corrected for beam, pixelization, and the additional 4' 
smoothing of the correlation function, and finally averaged 
into flat band-powers Ci with widths = 50 and posi- 
tions £ that of de Bernardis et al. (2000). Figure p| shows 
Ci for all realizations, together with the corresponding er- 
ror bars. They are to be compared with the optimal error 
bars computed from the approximate formula 



AC, = ^^^^^JC, + n.-Hii+l)W-V2n), (3) 

where fsky is fraction of the sky observed, is the noise 
variance and Wi is the beam of the experiment. According 
to Figure ^ ojir estimator is nearly optimal and unbiased. 
In any given f-band, the differences between the measured 
and theoretical mean and variance are at most 10% of the 
optimal variance. The MC method for determining the 
errorbars converges remarkably fast: a few tens of real- 
izations give an accurate estimate of both the mean and 
variance. For real data, the smoothed measured C^'s them- 
selves should be taken for generating MC realizations. 

In addition, we have estimated the cross-correlations be- 
tween bands, and found them to be consistent with zero, 
decreasing as ~ 'i/^/N^ with the number of realizations 
Nr- These experiments show that our method does not 
induce spurious correlations into the measurements. 

The right panel of Figure |^ displays distributions of sev- 
eral recovered band-powers, illustrating different regimes 



of signal or noise domination of the total variance. We 
found that the distribution in all bands is well fitted with 
an offset lognormal probability density (BJK98), although 
it is indistinguishable from a Gaussian in most bands. 

In summary, our method has nearly identical mean, vari- 
ance, cross-correlations, and distribution as one would ex- 
pect from standard maximum likelihood methods. The 
present implementation of the above pipeline takes about 
25 minutes of serial CPU on a small workstation, including 
artificial map generation, measurement of the correlation 
function, and estimating the C;'s individually, and finally 
summarizing them into band-powers. The required time is 
dominated by the measurement of the two-point correla- 
tion function (about 20 minutes); obtaining the Cj's from 
it is negligible (few seconds). 

These measurements illustrate that pairwise estimators 
with heuristic weights are feasible, their speed is far su- 
perior to maximum likelihood techniques, and they can 
be rendered nearly optimal. While details of the above 
recipe are subject to honing for future practical applica- 
tions, our present code is capable of analyzing MAP on a 
sub-supercomputer class CPU without explicit reliance on 
symmetries; a feat no other method could claim. 

4. DISCUSSION 

We presented a novel method to extract Cis from large 
CMB maps via two-point correlation functions. We have 
cast our determination of into 3 steps: fine-grained C{9) 
estimation, the pixel-related 4' Gaussian smoothing of it, 
followed by multiplication and Gauss-Legendre integra- 
tion. The resulting expression for Cg can be itself under- 
stood as a sum over pixel pairs, with a pair-weighting pro- 
portional to the Gaussian smoothing of P^. In this respect 
our estimator belongs to the class of quadratic estimators 
for Cf, of which other examples (e.g., BJK98) have been 
used in the past and are expected to give similar results. 

In the present example the measured errorbars were at 
most 10% larger then the theoretically minimal ones; this 
negligible suboptimality allowed the reduction of the CPU 
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time from months (Borrill 1999, Table 3) to hours. For 
future large missions the difference will be even more dra- 
matic: optimal methods would improve the errorbars only 
slightly at a prohibitive cost. 

Our new approach has several advantages in compar- 
ison to previous techniques. It scales as N"^ even in its 
most straightforward implementation, compared to typ- 
ical quadratic estimators which scale as . This scal- 
ing can be further improved up to N log N . While in the 
present implementation we used diagonal noise, no other 
(e.g., asymuthal) symmetries were assumed about noise, or 
geometry of the map. Cut out holes around bright sources, 
galactic cut, any irregularity in the sampling, make no dif- 
ference in the speed or performance of the algorithm. To 
illustrate this, Equation ^ was used to determine the cor- 
relation function from COBE DMR data, which has inho- 
mogeneous noise. The recovered Ci's are consistent with 
those obtained from implementations of quadratic estima- 
tors discussed in BJK98. In contrast with any other previ- 
ous attempt, our approach is straightforward to generalize 
for non-diagonal noise matrix (see below). 

While our method in its present form is already practical 
for analyzing megapixel CMB maps, it has the potential 
for further generalizations. General non-diagonal noise ap- 
pears to pose a serious problem. Fast iterative map making 
methods (Wright et al. 1996, Prunet et al. 2000) capable 
of handling large data sets only furnish the weight matrix, 
Wij — N~j^. Since inversion of the noise matrix is again 
an N'^ problem, our estimator of Equation might not be 
directly applicable. Instead, we propose a new estimator 
using MC realizations of artificial noise 

M 

C(cos0) = ^ /,,(A,A, - - ^ nfn^^), (4) 
I'j k=i 

where is one of M realizations of the noise for pixel 
k. This plays the role of a MC inversion of the weight 
matrix Wij = . Generation of in the time domain. 



where it has simple correlation structure, is straightfor- 
ward. Its iterative projection into pixel space is equivalent 
to map-making. This is feasible even when storing the 
noise-matrix would be prohibitive, as for Planck. 

To further improve the performance of our method, 
we will implement a minimum variance weighting scheme 
(Feldman Kaiser, & Peacock 1994, and Colombi, Szapudi, 
& Szalay 1998). This corresponds to down- weighting mea- 
surements with their variance, and might be important 
in maps, where various pairs contributing to the corre- 
lation function have widely differing errors. Otherwise, 
the uniform weighting scheme is nearly optimal (Colombi, 
Szapudi, & Szalay 1998). The present simple weighting 
scheme can be improved heuristically via individual pixel 
weighting reflecting differences in sky coverage. More com- 
plex pair weighting can be defined iteratively, using the 
MC estimates of the variances. 

The present Letter used a simple N"^ code to calculate 
the correlation function, and this is perfectly adequate for 
LDB surveys with N ~ 10^ pixels, and, with supercom- 
puters, even for MAP. Nevertheless, the development of 
an TV log code is under way (Szapudi & Colombi 2000, 
Connolly, Nichol, Moore, & Szapudi 2000). 

Our technique has reestablished the utility of corre- 
lation functions for CMB studies. This approach has 
further potential applications: it is naturally generaliz- 
ablc for the assessment of non-Gaussianity in CMB maps 
via fc-point correlation functions, with implementations as 
fast as iV(logA)'^~^ (Sunyaev-Zeldovich effect, lensing of 
CMB); it is equally useful for polarization correlation func- 
tions and for obtaining the appropriate CiS from them. 
We have found that the statistical information is con- 
densed into singular features of the correlation function, 
(see also Bashlinsky & Bertschinger, in prep.). This sug- 
gests that direct parameter estimation from the CMB two- 
point correlation function might be fruitful as well. 

Other applications include correlations of the infrared 
(SCUBA, BLAST, SIRTIF, FIRST) and optical back- 
ground, and weak gravitational lensing. 
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